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Abstract 

We study the three-dimensional static configurations of nonhomogeneous Kirch- 
hoff filaments with periodically varying Young's modulus. This type of variation may 
occur in long tandemly repeated sequences of DNA. We analyse the effects of the 
Young's modulus frequence and amplitude of oscillation in the stroboscopic maps, 
and in the regular (non chaotic) spatial configurations of the filaments. Our analysis 
shows that the tridimensional conformations of long filaments may depend critically 
on the Young's modulus frequence in case of resonance with other natural frequen- 
cies of the filament. As expected, far from resonance the shape of the solutions 
remain very close to that of the homogeneous case. In the case of biomolecules, 
it is well known that various other elements, besides sequence-dependent effects, 
combine to determine their conformation, like self-contact, salt concentration, ther- 
mal fluctuations, anisotropy and interaction with proteins. Our results show that 
sequence-dependent effects alone may have a significant influence on the shape of 
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these molecules, including DNA. This could, therefore, be a possible mechanical 
function of the "junk" sequences. 
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The study of tridimensional structures of filamentary objects is of great in- 
terest in several areas of knowledge, ranging from microscopic to macroscopic 
systems. Examples of macroscopic systems in Engineering are the stability of 
suboceanic cables [1,2] and installation of optical fibers [3,4]; in Biology, the 
shape of climbing plants [5] ; and, in Physics and Mathematics, the shape and 
dynamics of cracking whips [6]. The microscopic examples are in the area of 
the Structural Biology, as in the study of the structure of biomolecules [7,8] 
and bacterial fibers [9,10]. 

Filamentary systems are usually modeled as thin uniform rods. However, 
nonuniformities in the filament properties can affect significantly its tridimen- 
sional structure, so in this work we study certain resonant variations in the 
shape of helical rods induced by periodic variations in the stiffness of the rod. 
One of the motivations for this theoretical study is, on one hand, the obser- 
vation that the DNA stiffness is sequence-dependent [11]. On the other hand, 
it is also known [12,13] that a substantial fraction of all eukaryotic genomes 
are composed of tandemly repeated sequences of base-pairs. These repetitive 
DNAs are formed by nucleotide sequences of varying length and composition 
reaching up to 100 megabasepairs of length [13]. Usually they are regarded as 
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"selfish" or "junk" DNA [14] because they seem to have little or no functional 
role. In fact, some studies suggest that the behavior of repetitive sequences can 
be, in some cases, beneficial to the organism and, in others, harmful [12,13]. 
It could also be related to some form of cancer [15]. 

In this work we consider rods with small periodic variation of the Young's 
modulus, motivated by the existence of a large number of repetitive sequences 
of DNA. So the numerical calculations presented here have been performed 
using DNA parameters. Nevertheless, the qualitative results remain valid for 
general rods with periodic stiffness variation. We remark that ideal elastic rod 
models are not considered to give realistic solutions for the spatial structure 
of the DNA or other biomolecules [16]. Therefore, our results for the tridimen- 
sional configurations must be considered as general mechanical tendencies due 
to sequence-dependent variations of the Young's modulus, rather than exact 
solutions for the DNA structure. Realistic models for the DNA consider base- 
pair parameters, as the recent theory of sequence-dependent DNA elasticity 
proposed by Coleman, Olson and Swigon [17], where the elastic energy de- 
pends on a function of the six kinematical variables relating the relative orien- 
tation and displacement between successive base-pairs. Thermal fluctuations 
play also an important role in the structure of these molecules and, there- 
fore, statistical mechanical approaches are more appropriate to model their 
spatial configurations [18,19,20]. Nevertheless, in order to analyse the effect 
of periodic nonhomogeneities in equilibrium configurations, it is sufficient to 
take into account thermal fluctuations just roughly. This has been done by 
adopting an excess of 5% of linking number [21]. 

A table with the Young's modulus of all 32 trinucleotide units was recently 
obtained in [22] that can be used as a reference for the amplitude of the varia- 
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tion. We analyse the effects of the frequency of these periodic modulations in 
the Young's modulus, in terms of dynamical stroboscopic maps and directly 
in the tridimensional structure of the rod. We are particularly interested in 
understanding how the sequence-dependent mechanical properties of the fil- 
ament can cause its shape to deviate from the well known uniform solution, 
namely, the helix. 

Mielke and Holmes [23] demonstrated that the variation of the bending stiff- 
ness along the rod can cause spatially complex tridimensional shapes and sen- 
sitivity with respect to initial conditions. They described infinitely long rods 
as initial value problems (IVP) and used analytical techniques based on dy- 
namical systems theory to study some specific hyperbolic fixed points related 
to homoclinic orbits. 

It should be stressed that here we analyze a case where it is not possible to 
use the perturbative methods of dynamical systems theory, therefore we had 
to resort to numerical simulations. 

The Kirchhoff rod model has been extensively used in the literature to model 
continuous rods [24,25,26]. For example, Shi and Hearst [27] and Nizette and 
Goriely [28] obtained and classified, respectively, all the solutions of the static 
Kirchhoff equations for homogeneous rods with circular cross section. Cole- 
man et al [29] made a complete analysis of the stability of DNA within the 
framework of the Kirchhoff rod model, subjected to strong anchoring end con- 
ditions. Recently, da Fonseca and de Aguiar [30] compared homogeneous and 
nonhomogeneous rods, with varying stiffness, subjected to boundary condi- 
tions. The effects of nonhomogeneous mass distribution in the dynamics of 
unstable closed rods have been analyzed by Fonseca and de Aguiar [31]. Yang 
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et al [8] have studied a particular case of nonhomogeneous Young's modulus 
for closed rods. Manning et al [32] have incorporated into the Kirchhoff model 
the sequence-dependent discrete data of the intrinsic curvature of DNA. In the 
present study we assume that the DNA is intrinsically an untwisted straight 
rod, but we include sequence-dependent stiffness. 

In the Kirchhoff model for an inextensible rod the Young's modulus appears 
in the equations through its bending coefficient. As there is no table with the 
bending coefficients for all di or trinucleotides, we shall consider the data in 
[22] as a reference for our calculations, as mentioned. 

We shall consider the Hamiltonian formulation of the Kirchhoff's equations. 
We shall follow the simplest derivation by Nizette and Goriely [28] instead of 
the rigorous Hamiltonian formulation by Mielke and Holmes [23]. The reader 
is referred to [25,26,31] for a derivation of the Kirchhoff model, and to ref. [28] 
for the Hamiltonian formulation. Hamilton's equations for the Kirchhoff model 
are analogous to those of a symmetric spinning top in a gravity field, with the 
arc length s along the rod playing the role of time. The main advantage of a 
Hamiltonian formulation is that the theory of chaotic Hamiltonian systems and 
stroboscopic maps can be directly applied to understand the spatial behavior 
of the filament. The length of the tandemly repeated sequences can reach up 
to 100 megabasepairs [13] while the length of the repeats is no more than a 
few hundreds of basepairs long. So we shall solve the Kirchhoff's equations 
as an initial value problem (IVP) to find the conformational solutions of the 
filament problem. 

The Hamiltonian for an elastic rod with circular cross section, in Euler angles, 
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where E = E(s) and /x = /i(s) are the scaled Young's and shear moduli, 
respectively. F is the intensity of the total contact force (constant) exerted on 
the cross section at s. The units used here are the same of refs. [25,26,31]. T 
varies between 2/3 (incompressible material) and 1 (hyper-elastic material), 
and it should be remarked that T has no influence in the equilibrium solutions. 
The momenta are defined by 

Pe = E{s)6' , 

P^ =T Ql i(s)(<t/ + il)'coae) , ( 2 ) 
P^ = E{s)ip' sin 2 9 + cos 9 , 

where the prime indicates the derivative with respect to s. 

In the case of a homogeneous filament, E(s) = 1 and /i(s) = 1 in the equation 
(1), and the Hamiltonian is written as: 

H = Y + §- +V{d) > (3) 

where the potential V(9) is: 

m = ij^5g!f + Fmt . (4) 

If P^ = Pfa then 9 = (straight rod) is an equilibrium solution, and if F ^ 
there is a second equilibrium solution (9 ^ 0) corresponding to a helix [33]. 
If P^p = —P<f>, the equilibrium solution is 9 = ir, and if F ^ there is a 
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second equilibrium solution as above. If \P^\ ^ \P^\ then V(9) (4) has a single 
minimum corresponding to the well known helix solution. Denoting by 9q the 
point of minimum of the potential (4), the frequency iv Q of small oscillations 
around 9o is given by: 

u = P$ + 2V(e )-6Fcos9 . (5) 

In the Appendix we show that the potential V{9) (4) cannot be well approxi- 
mated by a second order expansion in the neigborhood of its minimum (it is 
expanded up to order 6). 

The Euler angles 9, and ip connect a fixed Cartesian basis {ei,e 2 ,e 3 } to 
the local orthonormal basis dj = dj(s, t), i = 1,2,3, attached to each point 
of the rod. The direction of e 3 is chosen to be the direction of the constant 
force F. d 3 is chosen to be tangent to the curve x(s) that defines the axis of 
the filament, and di and d 2 are in the direction of the principal moments of 
inertia of the cross section (perpendicular to d 3 ). The momentum P^ is the 
e 3 -component of the angular momentum with respect to the axis of the rod 
and Pff, is the torsional moment, i. e., the momentum with respect to d 3 [28]. 
They remain constant along the rod even if E and /i depend on the arc length 
s. The Hamiltonian, eq.(l), will depend on s through E(s) and fJ,(s). 

We shall consider the following periodic variation of the scaled Young's mod- 
ulus: 

E(s) = 1 + a cos lus . (6) 

The parameter a is the amplitude of the Young's modulus periodic variation, 
and iv is the frequency of the oscillation. We are concerned only with the 
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the tridimensional shape of the rod and it should be stressed that the shear 
modulus fj,(s) does not affect the tridimensional configuration. 

To obtain the equilibrium configurations we first solve Hamilton's equations 
for 9 and Pq. Then, we solve Eq.(2) for ip and reconstruct the filament by 
integrating d 3 along s: 

s 

x(s) = J [(sin 9 cos ip) ex + (sin 9 simp) e 2 + (cos 9) e 3 ]ds' (7) 
o 

x(s) is a function of the initial conditions 9(s = 0) = 9 and Pq(s — 0) = P . 
Without lack of generality, P can be set equal to so that 9 will be a 
conformation parameter. In solving the equation for ip we set its initial value 
^o = 0. 

In what follows we present numerical calculations performed with the following 
fixed mechanical parameters: a = 0.1, P^ = 0.086, P^ = 0.043 and F = 20pN. 
These parameters, excepting the force F, are written in properly scaled units. 
The maximum value for a is 0.66, in accordance with the table of the DNA 
Young's modulus presented in reference [22]. The value of P^ used corresponds 
to an excess of 5% of the linking number [21] due to thermal fluctuations. The 
value of the force corresponds to a compressing force consistent with the values 
in the literature [34]. 

Fig .1 displays nine stroboscopic maps on the 9 — Pq plane for different values 
of frequency u. We start with uo = 0.60cj (Fig .la) where a larger stability 
island encloses the main equilibrium point at 9 ~ 2.08rad and P e = 0, and 
a smaller island is seen on the left, at 9 ~ 0.5rad. The frequency goes up 
to uo = 2.00c<;o- We recall that spatial chaos has been observed before in the 
Kirchhoff equations [23,35]. 
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As we go through the sequence of stroboscopic maps displayed in the Fig .1, 
the large island in the plate (a) slowly shrinks and eventually disappears at 
uj ps uj (Figs .le and .If). A second important island appears at uj ps 0.82u 
(Fig .lb). This new island increases in size and moves towards the right as uj is 
increased. Besides these two main islands, a number of smaller and short-lived 
islands pop up and disappear as uj changes, a phenomenon typical of chaotic 
maps. We shall concentrate our study on the two main islands described above, 
since they dominate the stroboscopic maps and last for large intervals of uj. 

We shall now investigate the differences in the shape of the tridimensional 
configurations corresponding to the two equilibrium points lying at the center 
of these islands. In order to construct the rods we solved the Hamiltonian 
equations using the values of the equilibrium point for 6 and P e as initial 
conditions 6$ and P$ and used the equation (7) to construct the filament. 

The tridimensional configuration corresponding to the equilibrium point changes 
as the frequency is varied. The shape evolution is displayed in Figs .2 and .3 
for the two main equilibrium points mentioned above. 

In Fig .2, panels (a) to (d), we show the shape evolution of the configuration 
corresponding to the first main equilibrium point which lies in the center of 
the main island appearing in the Fig .la, and in the center of the island on 
the right in the Figs .lb, .lc and .le, respectively. We can see that the shape 
of the rod deviates more and more from the helix pattern as uj is increased, 
becoming rather twisted for uj = 0.92cj . 

In Fig .3, panels (a) to (d), we show the shape evolution of the configuration 
corresponding to the other main equilibrium point which lies in the center 
of the 'left island' (born in the Fig .lb). The configurations shown in the 
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Figs .3a-d correspond to the frequency values use in the Figs .lb, .lc, .If and 
.lh, respectively. The behavior of this sequence of rod shapes is the reverse of 
that corresponding to the first equilibrium point (Fig .2). As u; increases, the 
shape becomes less coiled and eventually recovers the near-helix shape, similar 
to the rod in Fig .2a, corresponding to the first equilibrium point (Fig .la). 

Finally, when uj = 2co , Fig .li, a 'period-doubling' bifurcations occurs. The 
orbit at the center of the island becomes unstable and a new stable equilibrium, 
with twice the original period, appears. 

The sensitivity of the shape of the nonhomogeneous rods to the amplitude of 
the nonhomogeneity can also be tested. Fig .4 shows, for the same mechanical 
parameters of the previous figures, the helix solution of the homogeneous case 
(left), the solution for a = 0.001 (middle) and the solution for a = 0.01 (right), 
in the resonant case u = u . Also, these solutions can be compared to that in 
the Fig .3c. We can see that even for very small values of the amplitude a, 
the tridimensional configuration deviates fast from the helix solution at the 
resonance. 

It is interesting to notice that, as the frequency increases, the position of the 
equilibrium points move in the direction of increasing 9. Fig .5 displays the 
value of 6 corresponding to the equilibrium point related to the main (circles) 
and to the left (square) islands as function of uj. The dotted line indicates 9 
which is the position of the equilibrium point of the potential V{6) related to 
the homogeneous case. As we can see in the Figs .2a and .3d, the shape of 
the corresponding tridimensional configuration becomes similar to the helix 
when the position of the equilibrium position gets close to 9 (homogeneous 
case equilibrium point position). 
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The main result of this numerical experiment is that the tridimensional confor- 
mations of long filaments may depend critically on sequence-dependent prop- 
erties if these are in resonance with other natural periods of the filament. As 
expected, in the limit of very low or very high frequencies, as compared to 
ujq, the shape of the solutions remain very close to that of the homogeneous 
case. In the case of biomolecules, it is well known that various other elements, 
besides sequence- dependent effects, combine to determine their conformation, 
like self-contact, salt concentration, thermal fluctuations, anisotropy and inter- 
action with proteins. Our results show that sequence-dependent effects alone 
may have a significant influence on the shape of these molecules, including 
DNA. This could, therefore, be a possible mechanical function of the "junk" 
sequences. 

This work was partially supported by the Brazilian agencies FAPESP, CNPq 
and FINER 

Appendix 

Here we show that the potential V(9) of equation (4) cannot be well approx- 
imated by an expansion up to order 2 or 3 around its minimum at (0 = O ). 
This is due to the presence of sin 2 in the denominator of one of the terms of 
the V(9). To illustrate this, we expand V(9), up to order 6 in (9 — 9 ), for the 
same numerical parameters used in this paper, 

V(9) = V(9 ) + 0.02 05(0 - O ) 2 + 0.0227(0 - O ) 3 + 0.0235(0 - O ) 4 + 

(•1) 

0.0257(0 - O ) 5 + 0.0275(0 - O ) 6 + 0[(9 - O ) 7 ] , 
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where 9 ^ 2.043 for this case. 

Since the coefficients of the terms (9 — 9 ) n have the same order of magnitude, 
it is necessary to check if (9 — 9 ) « 1 for all 9(s), i.e., along the rod. We 
found that for the frequency uj of the Young's modulus far from the resonance 
(\uj — uj \ » 0), the solutions corresponding to the equilibrium points have 
(9(s) — 9 )max — 0.03, where the subscript MAX means "maximum value 
for all s". But at the resonance, uj = u>o, the solution corresponding to the 
equilibrium point has (9(s) — #o)max — 0.5. (9(s) — 6> )max becomes even larger 
than 0.5 if we consider the solution related to the new equilibrium point (the 
new island that appeared in the map displayed in the Fig .lb). 

Therefore, the perturbative method of the dynamical systems theory is not 
applicable to analyzing this case. 
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Fig. .1. Stroboscopic maps for = 0.086, P<f, = 0.043 and F = 20pN. The frequency 
of the Young's modulus oscillation in each map is: (a) uj = 0.60cjo; (b) uj = 0.82wo; 
(c) u = 0.85w ; (d) w = 0.90w ; (e) w = 0.92w ; (f) w = l.OOwo; (g) w = 1.20wo; (h) 
c<j = 1.60wo; (i) oj = 2.00o;o. 
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Fig. .2. Shape evolution of the conformations corresponding to the center of the 
main island seen in Fig.. la. (a) u = O.GOujo (Fig.. la); (b) lo = 0.82c<jo (Fig.. lb); (c) 
co = 0.85w (Fig..lc); (d) u = 0.92u;o (Fig..le). 
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Fig. .3. Shape evolution of the conformations corresponding to the center of the 
left island that appears when u > 0.81^0 (Fig.. lb), (a) to = 0.82wo (Fig.. lb); (b) 
u = 0.85lj (Fig..lc); (c) u = w (Fig.. If); (d) uj = 1.60wq (Fig..lh). 
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Fig. .4. Tridimensional shape of the configurations corresponding to the center island 
that appears in the stroboscopic maps for to = ujq and different a (stroboscopic maps 
not shown). From left to right, homogeneous case, a = 0.001 and a = 0.01. 
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Fig. .5. Position of the equilibrium points as function of the frequency lo, circles 
corresponding to the center of the main island in the stroboscopic maps of Fig .la-e, 
and squares to the center of the left island that appears when to > 0.81u;o (Fig.lb-h). 
The dotted line corresponds to 9q, homogeneous case equilibrium point position. 
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